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Abstract We use a recently developed technique, which allows to perform few Lanczos 
steps on a given wavefunction even for large system sizes, to investigate the 
t—J model in the physical parameter region and to check the stability of the 
BCS d-wave variational wavefunction [1]. Our statistical Lanczos algorithm, 
which extends and improves the one Lanczos step proposed in Ref. [2], has 
been extensively tested on the small L = 18 sites cluster where many Lanczos 
iterations can be performed exactly. In this case, at doping 8 ~ 10% the 
BCS wavefunction represents a very good initial state to achieve extremely 
accurate energies and correlation functions with few Lanczos iterations. For 
large sizes (L < 98) the behavior is similar: the low-energy d-wave order 
parameter Pa is weakly affected by a couple of Lanczos iterations in the low 
doping 8 ~ 10% region, whereas the energy is considerably lowered. As a 
further test of our calculation we have computed the variance of the Hamiltonian 
AE P = ((H 2 ) - (H) 2 )/L 2 on the BCS wavefunction with p = 0,l,2Lanczos 
steps. For large p, when the Lanczos algorithm converges to the exact ground 
state, the variance vanishes exponentially with increasing p. The remarkable 
reduction of the variance, observed for p = 1,2 Lanczos steps even for the 
largest lattice size considered, suggests a smooth and rapid convergence to the 
exact ground state. These results support the existence of off-diagonal long-range 
d-wave superconducting order in the two-dimensional t—J model. 

1. INTRODUCTION 

One of the most important question raised after the discovery of high-Tc 
superconductivity is whether a simple model of strongly correlated electrons 
can capture the low-energy physics of real materials. In particular it is still a very 
much debated issue whether a purely repulsive electronic interaction can give 
rise to a d-wave superconducting ground state by doping an antiferromagnetic 
Mott insulator with a small amount of holes. The simplest model that has been 
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proposed immediately after the discovery of high-Tc superconductors is the 
t—J model [3,4] 

H = J • Sj - jfiih^j - t ^2 c\ a Cj, a , (1.1) 

where c| CT = c| CT (1 — n^), n« = ^2 a fii y(T is the electron density on site i, 

Sj = Sctct' 4o- r o-,o-'Ci,o-' is the spin operator and T a y are Pauli matrices. In 
the following we put t = 1. 

After many years of intense numerical and theoretical efforts there is no 
general consensus on the properties of this simple Hamiltonian and of the re- 
lated Hubbard model. From the numerical point of view, the density matrix 
renormalization group (DMRG) [5] predicts that a charge density wave insta- 
bility [6], nowadays called striped-phase for its one dimensional character, is 
strongly competing with pairing and superconductivity [7]. Within DMRG 
it appears therefore difficult to explain Copper-Oxide superconductors with 
a simple one-band model, especially because, in order to be consistent with 
photoemission experiments [8], a negative next-nearest-neighbor hopping am- 
plitude t' has to be included in the model. In fact, the negative t' suppresses 
even further superconductivity, so that the t—J model becomes unrealistic to 
describe the low-energy physics of high-Tc superconductors. However, the 
DMRG results, though quite accurate, are not exact in two dimensions. More- 
over, for technical reasons it is possible to consider only particular boundary 
conditions (open in one direction and periodic in the other), which certainly 
make the DMRG calculation still far to be representative of the thermodynamic 
limit. 

Quantum Monte Carlo (QMC) is an appealing alternative numerical ap- 
proach. This numerical method is still severely limited for two-dimensional 
fermionic systems by the well-known sign problem and is consequently biased 
by the initial guess of the ground state used to control this numerical instability 
[9, 10, 11]. However this technique has the important advantage to work very 
well with periodic boundary conditions since translation invariance can be ex- 
plicitly used to improve the efficiency of QMC algorithms. In particular, an 
approximate ground state can be obtained starting from a translation invariant 
wavefunction \ipc) by applying exactly few powers of the Hamiltonian (—H) p 
[12] or many approximate ones (by using for instance the Fixed-Node FN ap- 
proximation [9]). However, within various QMC schemes the situation is still 
controversial. First Heeb [2] and Khono [13] have found d-wave superconduc- 
tivity in a reasonable parameter range. Later Shih and co-workers [14], using a 
very similar method, have excluded drastically this possibility. The latter results 
were obtained within the monotonic-behavior assumption of the off-diagonal 
superconducting order parameter as a function of the number p of Hamiltonian 
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powers applied to the initial wavefunction. This assumption, although reason- 
able, is highly questionable. By contrast, in a recent QMC work [15], a clear 
tendency to d-wave superconductivity in the t—J model was found. Moreover, 
very recently, an almost realistic phase diagram with a corresponding high-Tc 
d-wave superconducting transition has been obtained for the Hubbard model 
within the Dynamical Mean Field Approximation [16, 17]. Furthermore, in 
the contest of the Hubbard model, also a weak-coupling renormalization group 
approach [18, 19] gives rise to a d-wave order parameter in a large region of 
the phase diagram. The latter results strongly support the relevance of a single 
band model for the explanation of high-Tc superconductivity. 

2. NUMERICAL METHOD 

In this work we make a further attempt to clarify the controversial numerical 
findings on the issue of d-wave superconductivity in the t—J model, using 
the statistical few Lanczos-step technique (FLST), efficiently implemented by 
means of the stochastic reconfiguration (SR) [10, 1 1]. Within the latter scheme 
all kind of correlation functions can be computed efficiently without any mixed 
average [11] approximation: an enormous advantage compared to the FN or 
to the original SR technique [10, 11]. In these cases this bias can be removed 
at the price of adding a small field coupled to the desired correlation function 
[15]. However, it turns out that it is extremely difficult and computationally 
demanding to work in the small field limit when unbiased correlation functions 
can be obtained. FLST is instead a very good compromise that solves efficiently 
this numerical problem of QMC methods. 

The wavefunctions that we are able to sample statistically read: 

Wp) = (i + E«^ p )|V>g> < l2 ) 
k=l 

with parameters for A; = 1, • • • ,p minimizing the energy expectation 

value (ip p \H\tfjp} / (t/jpli/jp) . For any p it is simple to show that the wavefunc- 
tion (1.2) corresponds exactly to apply p Lanczos step iterations to the initial 
wavefunction |V>g)- This wavefunction is sampled statistically with the SR 
technique, by using in the reconfiguration scheme the first p powers of the 
Hamiltonian. Unlike the previous method [10, 11] the reference wavefunction 
) is not evolved during the Markov iteration. | i/j* ) is instead kept statistically 
equal to the initial wavefunction \ipc) using the variational scheme proposed by 
Hellberg and Manousakis [20], which highly reduce the statistical fluctuations 
related to the SR technique. The equivalence of FLST to the standard Lanczos 
algorithm will be discussed in a forthcoming paper [21]. 



4 



The initial wavefunction to which FLST will be applied can be written as 
follows [1]: 

\^ G ) = ^ P=0 ) = P P N J\D). (1.3) 

where \D) is a BCS wavefunction, which is an exact eigenstate of the following 
Hamiltonian: 

Hbcs = Ho + ^f^(^ + A) (1.4) 

At = E M ^(<t4 + 4,T<|) d-5) 
M 

where Hq = e k cj, a Ck, a is the free electron tight binding nearest-neighbor 

Hamiltonian, e& = — 2(cosfe E + cosk y ) — \i, \x is the free-electron chemical 
potential and A* creates all possible nearest-neighbor singlet bonds with d- 
wave symmetry being Mjj +1 or —1 if the bond (i, j) is in the x or y direction, 
respectively. Pjy and Pq are the projectors over the subspaces with a fixed 
number N of particles and no doubly occupied states. Finally the Jastrow factor 
J = exp ^7/2 ~~ j^ifij^ couples the holes via the density operators 

hi and contains another variational parameter 7 ~ 1 which scales an exact 
analytic form, obtained by approximating the holes with hard-core bosons at 
the same density, and applying the spin-wave theory to the corresponding XY 
model [22]. We note here that by performing a particle-hole transformation on 
the spin down c\ , — > (— l) l Cj,j., the ground state of the BCS Hamiltonian is 
just a Slater-determinant with N = L particles [23]. This is the reason why 
this variational wavefunction can be considered of the generic Jastrow-Slater 
form, a standard variational wavefunction used in QMC. Using the particle-hole 
transformation, it is also possible to control exactly the spurious finite system 
divergences related to the nodes of the d-wave order parameter. 

3. NUMERICAL TESTS 

In this section we show the accuracy of FLST applied to the BCS wavefunc- 
tion (1.3) on a small 18-site cluster, where exact results are available. 
Our main task is to compute the order parameter at finite system size 

P, = i« +2 |At|<), (1.6) 

where \ ipp ) and \ipp +2 } are the states with N and N + 2 particles, respectively. 
If Pd is finite in the thermodynamic limit this necessarily implies off-diagonal 
long-range order in the ground state. Following Ref. [15], it is convenient with 
an approximate technique to calculate a short-range quantity like P^, instead 
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Figure 1.1 Ground-state energy per site Eq/L and d-wave order parameter Pd as a function 
of the variance for iV = 16, L = 18, J = 0.4: full dots (optimal Abcs), empty dots 
(Abcs ~ > 0). Dashed lines are quadratic fits of the estimates with p = 0, 1, 2. 

of the more conventional long-range expectation value (ip p |AA^|V>^)/L 2 . In 
Table 1.1 we show a comparison between FLST and the exact results for 18 
and 16 electrons at J = 0.4. 

In Table 1 .2 we show Pd as a function of the number of Lanczos step iterations 
for the 18-site cluster at J = 0.4. In the same Table we have computed also 
the variance AE P = ((ifj p \H 2 \ip p } — {ij)p\H\il)p) 2 ) / L 2 , the overlap squared 
Zp = KV'plV'o)) 2 °f the FLST wavefunction with the true ground state \ipo), 
and the average sign of the FLST wavefunction: 

(S p ) = ^(x|^) 2 Sgn((x|^)(x|^o)) , (1.7) 

X 

where \x) denotes configurations with definite electron positions and spins. 
For an exact calculation, namely p » 1, both Z p — > 1 and S p — > 1, whereas 



N Abcs P S /L(FLST) Pd (FLST) E /L (Exact) P d (Exact) 



18 


0.55 


1 


-0.4765(1) 




-0.47668 




18 


0.55 


2 


-0.4775(1) 




-0.47749 




16 


0.20 


1 


-0.6541(1) 


0.1074(4) 


-0.65420 


0.10730 


16 


0.55 


2 


-0.6583(1) 


0.122(1) 


-0.65826 


0.12135 



Table 1.1 Comparison between the estimates of the ground-state energy per site Eo/L and of 
the d-wave order parameter Pd obtained with the exact and the statistical (FLST) application of 
p Lanczos steps on the variational wavefunction of Eq. (1.3). L = 18, N = 16, 18 and J = 0.4. 



AE p — > 0. The variance thus represents a very important tool to estimate the 
'distance' from the exact ground state when the latter one is not known. In 
particular whenever Z p ~ 1 the energy approaches the exact result linearly with 
the variance AE p , allowing us to estimate the error in the variational energy. 

This can be achieved by plotting the variational energies E v as a function 
of the corresponding variance AE P , and performing a very stable linear or 
quadratic fit to the AE p = exact limit (see Fig. 1.1). Similar fits can be 
attempted for correlation functions though, in this case, also a term oc ^/AE P 
is expected for AE P — ► 0. This term is however negligible for quantities like 
Pd that are averaged bulk correlation functions in a large system size (see the 
Appendix). In practice even in the small 18-site cluster the non-linear term 
turns out to be negligible (see Fig. 1.1). We believe that, being the convergence 
of the Lanczos algorithm particularly well behaved and certainly unbiased, the 
variance extrapolation method is in this case particularly useful and reliable. 
However for bad initial wavefunction (e.g., randomly generated) or very large 
sizes the approach to zero of the variance may behave rather wildly, requiring 
many Lanczos steps to reach the regime where the extrapolation is possible. 

As shown in Table 1.2 the quality of the variational BCS wavefunction (1.2) 
is exceptionally good, especially in the doped N = 16 case. Here Z p is larger 
than 0.9 even at the simplest p = variational level, and is drastically improved 
with really few Lanczos step iterations. Remarkable is also the behavior of the 
average sign S p which measures directly the accuracy of the BCS wavefunction 
phases, without caring about the amplitudes. In the undoped case the signs of 
the BCS wavefunction (So) can be proven to be exact, i.e., (So) = 1, having the 
BCS state the well-known Marshall signs, i.e., the phases of the exact ground 
state of Heisenberg model. For the two-hole case, the BCS nodes change in 
a non trivial way. Nevertheless, (So) remains very close to 1 and it is much 
higher than the average sign of the corresponding Gutzwiller wavefunction 
(Abcs —> 0), also shown in the table for comparison. 
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N 


&-BCS 


P 


(S P ) 


Zp 


AE p x I? 


E /L 




18 


0.00 





1.0000 


0.6898 


1.194 


-0.43833 


16 


0.00 





0.9656 


0.8306 


1.174 


-0.63847 


18 


0.80 





1.0000 


0.8850 


0.335 


-0.46639 


18 


0.80 


1 


1.0000 


0.9915 


0.042 


-0.47662 


18 


0.80 


2 


1.0000 


0.9995 


0.004 


-0.47752 


18 


0.80 


3 


1.0000 


1.0000 


0.0003 


-0.47759 


18 


0.80 


4 


1.0000 


1.0000 


0.00002 


-0.47759 


18 


0.80 


oo 


1.0000 


1.0000 


0.0 


-0.47759 


16 


0.55 





0.9891 


0.9260 


0.940 


-0.64559 


16 


0.55 


1 


0.9988 


0.9814 


0.191 


-0.65638 


16 


0.55 


2 


0.9999 


0.9942 


0.060 


-0.65826 


16 


0.55 


3 


1.0000 


0.9983 


0.018 


-0.65882 


16 


0.55 


4 


1.0000 


0.9995 


0.005 


-0.65898 


16 


0.55 


5 


1.0000 


0.9999 


0.002 


-0.65902 


16 


0.55 


6 


1.0000 


0.9999 


0.0005 


-0.65904 


16 


0.55 


7 


1.0000 


1.0000 


0.0001 


-0.65904 


16 


0.55 


oo 


1.0000 


1.0000 


0.0 


-0.65904 



Table 1.2 Average sign {S p ), overlap squared on the exact ground state Z p and variance times 
the volume squared AE P x L 2 obtained applying exactly p Lanczos steps on the variational 
wavefunction of Eq. (1.3). L = 18, TV = 16, 18 and J = 0.4. 
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These results suggest that there is a tendency to d-wave BCS pairing in the 
t—J model at ~ 10% doping and J ~ 0.4, and that the BCS wavefunction is a 
particularly accurate wavefunction to describe the small and even zero doping 
region of the t—J model. 

4. LARGER SIZE CALCULATIONS 

Though few Lanczos steps may appear inadequate for large system size, 
this simple scheme is instead providing us very good variational energies up 
to L ~ 100 sites in the t—J model, even when this variational energies are 
compared with more complicated schemes like the FN. In order to show that 
the Lanczos scheme remains effective for larger sizes it is useful to consider 
first a relevant case where a numerically exact solution is possible: the zero 
doping limit of the t—J model, i.e., the Heisenberg model. In Fig. 1.2 we 
plot the energy results calculated for few Lanczos steps as a function of the 
variance. The fact that the energy is approaching the exact result smoothly 
with the variance both in the 50- and 98-site clusters (with a slightly larger 
curvature in the latter case) indicates that the projected BCS wavefunction 
should have a substantially large overlap squared Z p with the exact ground 
state of the Heisenberg model even on these large system sizes. This is also 
confirmed by the behavior of the square order parameter (Fig. 1.3), which 
is considerably improved by few Lanczos steps. This result is particularly 
important since the starting BCS wavefunction has not antiferromagnetic long- 
range order, whereas the two-dimensional Heisenberg model is widely believed 
to be antiferromagnetically ordered. 

The above test represents also a further strong evidence that the ground- 
state wavefunction of the Heisenberg model is smoothly connected to a d- 
wave BCS superconducting wavefunction. This circumstance represents a very 
interesting numerical fact that clearly supports the experimental observation 
of high-Tc d-wave superconductivity coming just upon a small doping of a 
quantum antiferromagnet. 

Within the Lanczos approximate states \ip p ) (p = 0, 1, 2) acting on the best 
variational BCS wavefunction, we have performed a finite-size scaling of the 
order parameter as defined in Eq. (1.6), at a fixed doping S = 13.3% 
(corresponding to N = 84 on the largest size L = 98). Since Pd is computed 
between two states with N and N + 2 particle numbers, we assume that it refers 
to the intermediate doping 5 = 1 — (N + l)/L. For smaller sizes the doping 
8 = 13.3% is not possible and we have interpolated linearly Pd between the 
two fillings closest to this doping. The main results of this paper is then shown 
in the Fig. 1.4. Here the size scaling for p = 0, 1, 2 Lanczos steps and for their 
zero-variance extrapolation clearly indicates a finite Pd in the thermodynamic 
limit. The stability of the BCS variational wavefunction is evident from this 
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Figure 1.2 Energy per site of the finite-size Heisenberg model. Comparison of exact results 
(indicated by arrows) and the approximate p = 0, 1, 2 Lanczos step iterations over the projected 
d-wave wavefunction. Continuous lines are quadratic fit of the data. 



figure: remains much larger than the corresponding value of the Gutzwiller 
(Abcs - > 0) metallic state. Instead, as shown in Fig. 1.5 and similarly to the 
undoped case, the energy is very much lowered by FLST, suggesting that the 
method remains effective even for large sizes and finite doping. 

5. CONCLUSIONS 

In conclusion we have provided evidences in favor of the stability of d-wave 
pairing correlations in the t—J model. Our result disagrees with a recent 
QMC calculation [14], for which, however, the self-consistent assumption 
mentioned in the introduction leads to a very poor, high-energy variational 
state. Furthermore, also the DMRG results provide a qualitatively different 
scenario, namely that stripes suppress d-wave superconductivity [7]. In order 
to clarify this disagreement, we have performed a simulation at J = 0.4 for a 
cluster with L = 6 x 12, 8 holes, and open boundary conditions on the long 
direction, as usually done within DMRG. Our accuracy is in this case worse 
than with periodic boundary condition, especially for quantities that are not 
bulk average correlation functions (see the Appendix). For P& instead we have 
found that, even at the simplest variational level open boundary conditions 
strongly suppress Pj by more than 30%, much more than the difference found 
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Figure 1.3 Order parameter m = v/ S(ty, tt) / L in the finite-size Heisenberg model (S(tt, it) 
being the spin isotropic antiferromagnetic structure factor). Comparison of exact results (in- 
dicated by arrows) and the approximate p = 0, 1, 2 Lanczos step iterations over the projected 
d-wave wavefunction. Continuous lines are quadratic fit of the data. Inset: finite-size scaling 
with the variational (BCS d-wave) wavefunction and with the variance extrapolated one. 



in Fig. (1.4). Thus we find that boundary conditions play a very important role 
and may strongly destabilize the uniform BCS d-wave ground state, as obtained 
within DMRG. This is consistent with experimental findings [24], showing 
that any small asymmetry in the CuO planes, acting in our point of view as 
a distorted boundary condition, may enhance the tendency to stripe formation 
and then suppressing superconductivity. Our results instead support the recent 
numerical works [16, 17, 18, 19], indicating that d-wave superconductivity can 
be obtained in a one band model with repulsive electron interaction. 

Appendix: Variance estimate of the error on bulk correlation 
functions 

In this Appendix we estimate the error on correlation functions assuming 
that the ground state \ipo) lS approximated with the wavefunction 

IV-o) = \i> P ) + e P W) (1-A.l) 

where (ipp\ip p ) = {ip'ty') = 1. an d W) represents a normalized wavefunction 
orthogonal to the exact one, (ipo\ip') = 0. We restrict our analysis to thermo- 
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Figure 1.4 Superconducting d-wave order parameter Pd in the t— J model as defined in the 
text in Eq. (1.6). These results were obtained starting by the variational wavefunction (triangles) 
defined in Eq. (1.3) and by applying to it one (squares) or two (dots) Lanczos steps. Black dots 
are obtained by quadratic extrapolation to the zero variance exact limit in order to estimate the 
size dependent error of the approximate variational calculations. The stars refer to the Gutzwiller 
wavefunction with Abcs — * and 7 = 1 in Eq. (1.3). Continuous lines are quadratic fit of the 
data. 



dynamically averaged correlation functions O, the ones which can be written 
as a bulk average of local operators Oc O = Y^iOi/L. This class of operators 
includes for instance the average kinetic or potential energy or the spin-spin 
correlation function at a given distance d, namely Oj = Sj • Sj + <2. If we use 
periodic boundary conditions the expectation value of Oj on a state with given 
momentum does not depend on i and the bulk average does not represents an 
approximation 

(MQ^M _ (MQbM _ c (1A2) 

(V'olV'o) (ipol^o) 

We show here that the expectation value of bulk-averaged operators O on the 
approximate state \ip p ) satisfy the following relation: 

(i> v \6\i> p ) = C + 0(e 2 p , ep/VZ), (1.A.3) 

thus implying that for large enough size the expectation value (1.A.3) ap- 
proaches the exact correlation function C linearly with the variance. The 
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Figure 1.5 Energy per site of the t — J model for 44 electrons on 50 sites (triangles) and 86 
electrons on 98 sites (squares), J = 0.4, for the approximate p — 0, 1, 2 Lanczos step iterations 
over the projected d-wave wavefunction Eq. (1.3). Continuous lines are quadratic fit of the data. 



validity of the above statement is very simple to show under very general 
grounds. In fact by definition: 



(if) p \d\i/) p ) =C + 2e p (i;'\d\^o) + 4(iP'\d\ip'). 



(1.A.4) 



The term proportional to e p in the above equation can be easily bounded by use 
of the Schwartz inequality: 



|(v/|o|Vo>| 2 = W\d - c|^ )| 2 < (Mo - c) 2 |Vo) 



(1.A.5) 



The final term in the latter inequality can be estimated under the general assump- 
tion that correlation functions C(d) = (?/>o|(Oi — C){Oi + d — C)|V'o)/(' i /'o|V ; o) 
decay sufficiently fast with the distance as a consequence of the cluster 
property: 



This concludes the proof of the statement of this Appendix, provided J2d C(^) 
is finite for L — > oo. 
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